Check correlation between Mn responsive and microbe responsive genes in ime course
library(tidyverse)
library(readxl)
library(edgeR)
get the Mn data
Mn <- read_csv("../output/Mn_deficiency.DEGs.all.anno.csv") %>%
select(AGI, logFC, FDR) %>%
filter(FDR < 0.05)
Missing column names filled in: 'X1' [1]Parsed with column specification:
cols(
X1 = [32mcol_double()[39m,
AGI = [31mcol_character()[39m,
logFC = [32mcol_double()[39m,
logCPM = [32mcol_double()[39m,
LR = [32mcol_double()[39m,
PValue = [32mcol_double()[39m,
FDR = [32mcol_double()[39m,
symbol2 = [31mcol_character()[39m,
full_name = [31mcol_character()[39m
)
Mn
Now get the microbe response data
load("../output/TimeCourseDGE_LRT.Rdata")
names(models)
[1] "day_time" "day" "time" "dge" "lrt" "DEG0.05"
DEgenes <- models %>%
mutate(toptags = map(lrt, ~ topTags(.,n=Inf)), #pull the DE gene info out of each model
toptags = map(toptags, ~ .$table), # get the table of DE info
toptags = map(toptags, ~ rownames_to_column(., "transcript_ID"))) %>%
unnest(toptags) %>%
group_by(transcript_ID) %>%
arrange(transcript_ID, FDR) %>%
select(day_time, transcript_ID, logFC, FDR)
DEgenes
Add Arabidopdis ID to each Brassica gene
First, get annotation
annotation <- read_csv("../../Brapa_microbes/Annotation/output/v3.0annotation/Brapa_V3.0_annotated.csv")
Missing column names filled in: 'X1' [1]Parsed with column specification:
cols(
.default = col_double(),
name = [31mcol_character()[39m,
chrom = [31mcol_character()[39m,
subject = [31mcol_character()[39m,
AGI = [31mcol_character()[39m,
At_symbol = [31mcol_character()[39m,
At_full_name = [31mcol_character()[39m,
At_gene_model_type = [31mcol_character()[39m,
At_short_description = [31mcol_character()[39m,
At_Curator_summary = [31mcol_character()[39m,
At_Computational_description = [31mcol_character()[39m
)
See spec(...) for full column specifications.
annotation
now join the annotation
DEgenes <- annotation %>%
select(name, AGI) %>%
right_join(DEgenes, by=c("name" = "transcript_ID"))
DEgenes
now add in the microbe data
gene_compare <- DEgenes %>%
inner_join(Mn, by = c("AGI"), suffix = c(".microbe", ".Mn"))
gene_compare
gene_compare <- gene_compare %>%
group_by(day_time) %>%
nest()
gene_compare
gene_compare <- gene_compare %>%
mutate(cortest = map(data, ~ cor.test(.$logFC.microbe, .$logFC.Mn, method="spearman", exact=FALSE)))
gene_compare %>% mutate(glance = map(cortest, broom::glance)) %>% unnest(glance) %>%
arrange(p.value) %>% select(day_time, estimate, p.value)
Plot the most significant one
gene_compare %>%
filter(day_time=="02-2_afternoon") %>%
unnest(data) %>%
ggplot(aes(x=logFC.Mn, y=logFC.microbe)) +
geom_point(alpha=.2)